Nature of phase transition in magnetic thin films 
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We study the critical behavior of magnetic thin hims as a function of the him thickness. We 
use the ferromagnetic Ising model with the high-resolution multiple histogram Monte Carlo (MC) 
simulation. We show that though the 2D behavior remains dominant at small thicknesses, there is 
a systematic continuous deviation of the critical exponents from their 2D values. We observe that 
in the same range of varying thickness the deviation of the exponent v is very small from its 2D 
value, while exponent f3 suffers a larger deviation. Moreover, as long as the film thickness is fixed, 
i. e. no finite size scaling is done in the z direction perpendicular to the film, the 3D values of the 
critical exponents cannot be attained even with very large (but fixed) thickness. The crossover to 3D 
universality class cannot therefore take place without finite size scaling applied in the z direction, 
in the limit of numerically accessible thicknesses. From values of exponent a obtained by MC, 
we estimate the effective dimension of the system. We conclude that with regard to the critical 
behavior, thin films behave as systems with effective dimension between 2 and 3. 

PACS numbers: 75.70.Rf Surface magnetism ; 75.40.Mg Numerical simulation studies ; 64.60.Fr Equilibrium 
properties near critical points, critical exponents 



I. INTRODUCTION 

During the last 30 years, physics of surfaces and ob- 
jects of nanometric size have attracted an immense inter- 
est. This is due to important applications in industry!^ 
An example is the so-called giant magneto-resistance 
(GMR) used in data storage devices, magnetic sensors, 
„i ?i 4 i?it In parallel to these experimental developments, 
much theoretical effort 7 -^ has also been devoted to the 
search of physical mechanisms lying behind new proper- 
ties found in nanoscale objects such as ultrathin films, 
ultrafme particles, quantum dots, spintronic devices etc. 
This effort aimed not only at providing explanations for 
experimental observations but also at predicting new ef- 
fects for future experiments. 

The physics of two-dimensional (2D) systems is very 
exciting. Some of those 2D systems can be exactly solved: 
one famous example is the Ising model on the square 
lattice which has been solved by Onsager^ This model 
shows a phase transition at a finite temperature T c given 
by sinh 2 (2 J/fcsT c ) = I where J is the nearest-neighbor 
(NN) interaction. Another interesting result is the ab- 
sence of long-range ordering at finite temperatures for 
the continuous spin models (XY and Heisenberg models) 
in 2D*£- In general, three-dimensional (3D) systems for 
any spin models cannot be unfortunately solved. How- 
ever, several methods in the theory of phase transitions 
and critical phenomena can be used to calculate the crit- 
ical behaviors of these systems?^ 

This paper deals with systems between 2D and 
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3D. Many theoretical studies have been devoted to 
thermodynamic properties of thin films, magnetic 
multilayers, ..i 7 ' 8 ' 12 ! 13 ! 14 In spite of this, several points are 
still not yet understood. We study here the critical be- 
havior of thin films with a finite thickness. It is known a 
long time ago that the presence of a surface in magnetic 
materials can give rise to surface spin-waves which are 
localized in the vicinity of the surfaced These localized 
modes may be acoustic with a low-lying energy or optical 
with a high energy, in the spin-wave spectrum. Low-lying 
energy modes contribute to reduce in general surface 
magnetization at finite temperatures. One of the conse- 
quences is the surface disordering which may occur at a 
temperature lower than that for interior magnetization. 16 
The existence of low-lying surface modes depends on the 
lattice structure, the surface orientation, the surface pa- 
rameters, surface conditions (impurities, roughness, ...) 
etc. There are two interesting cases: in the first case a 
surface transition occurs at a temperature distinct from 
that of the interior spins and in the second case the sur- 
face transition coincides with the interior one, i. e. exis- 
tence of a single transition. Theory of critical phenomena 
at surfaces 7 * 8 - and Monte Carlo (MC) simulation o 17 ' 18 of 
critical behavior of the surface-layer magnetization at the 
extraordinary transition in the three-dimensional Ising 
model have been carried out. These works suggested sev- 
eral scenarios in which the nature of the surface transition 
and the transition in thin films depends on many factors 
in particular on the symmetry of the Hamiltonian and on 
surface parameters. 

The aim of this paper is to investigate the effect of 
the film thickness on the critical behavior of the sys- 
tem. We would like to see in particular how the thick- 
ness affects the values of critical exponents. To carry 
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out these purposes, we shall use MC simulations with 
highly accurate multiple histogram technique ! 19 ! 20 ' 21 We 
consider here the case of a simple cubic film with Ising 
model. For our purpose, we suppose all interactions are 
the same even that at the surface. This case is the sim- 
plest case where there is no surface-localized spin-wave 
modes and there is only a single phase transition at a 
temperature for the whole system (no separate surface 
phase transition) 1 15 i 16 Other complicated cases will be 
left for future investigations. However, some preliminary 
discussions on this point for complicated surfaces have 
been reported in some of our previous papers. 22 ' 23 

The paper is organized as follows. Section II is devoted 
to a description of the model and method. Results are 
shown and discussed in section III. Concluding remarks 
are given in section IV. 



II. MODEL AND METHOD 
A. Model 

Let us consider the Ising spin model on a film made 
from a ferromagnetic simple cubic lattice. The size of 
the film is L x L x N z . We apply the periodic boundary 
conditions (PBC) in the xy planes to simulate an infinite 
xy dimension. The z direction is limited by the film 
thickness N z . If N z = 1 then one has a 2D square lattice. 

The Hamiltonian is given by 



H 



(id) 



Ui ■ (Tj 



(1) 



where 0$ is the Ising spin of magnitude 1 occupying the 
lattice site i, Y^Uj) indicates the sum over the NN spin 
pairs Oi and Oj. 

In the following, the interaction between two NN sur- 
face spins is denoted by J a , while all other interactions 
are supposed to be ferromagnetic and all equal to J = 1 
for simplicity. Let us note in passing that in the semi- 
infinite crystal the surface phase transition occurs at the 
bulk transition temperature when J s ~ 1.52 J. This point 
is called " extraordinary phase transition" which is char- 
acterized by some particular critical exponents. 17 ' 18 In 
the case of thin films, i. e. N z is finite, it has been the- 
oretically shown that when J s — 1 the bulk behavior is 
observed when the thickness becomes larger than a few 
dozens of atomic layers^ surface effects are insignificant 
on thermodynamic properties such as the value of the 
critical temperature, the mean value of magnetization at 
a given T, ... When J s is smaller than J, surface magne- 
tization is destroyed at a temperature lower than that for 
bulk spins' However, it should be stressed that, except 
at the so-called " extraordinary phase transition" ) 17 ' 18 the 
surface criticality has not been studied as a function of 
the film thickness. 



B. Multiple histogram technique 

The multiple histogram technique is known to repro- 
duce with very high accuracy the critical exponents of 
second order phase transitions ] 19 ' 20 ' 21 

The overall probability distribution^ at temperature 
T obtained from n independent simulations, each with 
Nj configurations, is given by 



P(E,T) 



Eg=i_gjCg) exp[£/fc B T] 



where 



eicp[f i }=J2P(E,T i ). 



(2) 



(3) 



The thermal average of a physical quantity A is then 
calculated by 



in which 



(A(T)) = '£AP(E,T)/z(T), 



(4) 



(5) 



Thermal averages of physical quantities are thus cal- 
culated as continuous functions of T, now the results 
should be valid over a much wider range of tempera- 
ture than for any single histogram. The xy linear sizes 
L = 20,25,30, ...,80 have been used in our simulations. 
We have tested that all exponents do not change in the 
finite size scaling with L > 30. So most of results are 
shown for L > 30 except for v where the lowest sizes 
L = 20, 25 can be used without modifying its value. 

In practice, we use first the standard MC simula- 
tions to localize for each size the transition tempera- 
tures Tq{L) for specific heat and T™{L) for suscepti- 
bility. The equilibrating time is from 200000 to 400000 
MC steps/spin and the averaging time is from 500000 
to 1000000 MC steps/spin. Next, we make histograms 
at 8 different temperatures Tj{L) around the transition 
temperatures T c f' m (L) with 2 millions MC steps/spin, 
after discarding 1 millions MC steps/spin for equilibrat- 
ing. Finally, we make again histograms at 8 different 
temperatures around the new transition temperatures 
Tq - m {L) with 2 x 10 6 and 4 x 10 6 MC steps/spin for 
equilibrating and averaging time, respectively. Such an 
iteration procedure gives extremely good results for sys- 
tems studied so far. Errors shown in the following have 
been estimated using statistical errors, which are very 
small thanks to our multiple histogram procedure, and 
fitting errors given by fitting software. 

We have calculated the averaged order parameter (M) 
(M: magnetization of the film), averaged total energy 
(E), specific heat C v , susceptibility \i first order cumu- 
lant of the energy Cjj, and n th order cumulant of the 
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order parameter V n for n 
are defined as 



1 and 2. These quantities 
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('•- = 
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(6) 
(7) 

(8) 

(9) 

(10) 



Plotting these quantities as functions of T for system 
size (L,N Z ), we can identify the transition temperature 
by looking at their respective behavior (maxima of C v 
and x, ...). Note that the transition temperatures for 
these quantities coincide only at infinite L. For large 
values of L, the following scaling relations are expected 
(see details in Ref. l2lh : 



U^ocL 1 ^, U 2 max cx L x l\ 

C™ ax = Co + C x L a l v 



and 



x max « 



(11) 



(12) 



(13) 



at their respective 'transition' temperatures T C (L), and 



Cu = Cu[T c (oo)]+AL 



-<xjv 



T c (oo) CX 



and 



T C (L) = T c (oo) + CaL-V", 



(14) 



(15) 



(16) 



where A, Co, Cl and C4 are constants. We estimate 
v independently from V7 nax and V 2 max . With this value 
we calculate 7 from ^ max and a from C™ ax . Note that 
we can estimate T c (oo) by using the last expression. Us- 
ing T c (oo), we can calculate [3 from Mxi^y The Rush- 
brooke scaling law a + 2(3 + 7 = 2 is then in principle 
verified. Finally, using the hyperscaling relationship, we 
can estimate the "effective" dimension of thin films by 
d e g = (2 — a) jv and the exponent 77 from the scaling law 
7 = (2 - rf)v. 

We note however that only v is directly calculated from 
MC data. Exponent 7 obtained from x max and v suffers 
little errors (systematic errors and errors from v) . Other 
exponents are obtained by MC data and several-step fit- 
ting. For example, to obtain a we have to fit C" lax of 



Eq. [12] by choosing Co, C\ and by using the value of v. 
So in practice, in most cases, one calculates a or /3 from 
MC data and uses the Rushbrooke scaling law to calcu- 
late the remaining exponent. However, for our precise 
purpose we will show in the following the results of all 
exponents v, 7, a and (3 calculated from MC data. We 
will show that the Rushbrooke scaling law is very well 
verified. The exponent a will allow us to estimate the 
effective dimension of the system. 



III. RESULTS 

We show now the results obtained by MC simulations 
with the Hamiltonian |T]). 

Let us show in Fig. [T] the layer magnetizations and 
their corresponding susceptibilities of the first three lay- 
ers, in the case where J s = 1. It is interesting to note 
that the surface layer is smaller that the interior lay- 
ers, as it has been shown theoretically by the Green's 
function method a long time agOi 15 ' 16 The surface spins 
have smaller local field due to the lack of neighbors, so 
thermal fluctuations will reduce their magnetization with 
respect to the interior layers. The susceptibilities have 
their peaks at the same temperature, indicating a single 
transition. 




FIG. 1: Layer magnetizations (a) and layer susceptibilities 
(b) versus T with N z = 5. 

Figure [2] shows total magnetization of the film and the 
total susceptibility. This indicates clearly that there is 
only one peak as said above. 
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FIG. 2: Total magnetization (a) and total susceptibility (b) 
versus T with N z = 5. 
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FIG. 3: (a) Susceptibility and (b) Vi, as functions of T for 
several L with N z = 11, obtained by multiple histogram tech- 
nique. 



Let us show now an example of excellent results ob- 
tained from multiple histograms described above. Figure 
[3] shows the susceptibility and the first derivative V\ ver- 
sus T around their maxima for several sizes. 

We show in Fig. 3] the maximum of the first derivative 
of In M with respect to (3 — (/cbT) -1 versus L in the 
In — In scale for several film thicknesses up to N z = 13. 
From the slopes of these remarkably straight lines, we 
obtain 1/v. We plot in Fig. [5] v as a function of thick- 
ness N z . We observe here a small but systematic de- 
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FIG. 4: Maximum of the first derivative of In M versus L in 
the In — In scale. 
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FIG. 5: v versus N z 



viation of v from its 2D value (v2D — 1) with increas- 
ing thickness. To show the precision of our method, we 
give here the results of N z = 1. For N z = 1, we have 
l/v= 1.0010 ± 0.0028 which yields v = 0.9990 ± 0.0031 
and j/u = 1.7537 ± 0.0034 (see Figs. M and below) 
yielding 7 = 1.7520 ± 0.0062. These results are in ex- 
cellent agreement with the exact results 1/2D — 1 and 
72/j = 1.75. The very high precision of our method is 
thus verified in the range of the system sizes used in the 
present work. 

We show in Fig. [S] the maximum of the susceptibility 
versus L in the In — In scale for film thicknesses up to 
N z = 13. We have used only results of L > 30. Includ- 
ing L — 20 and 25 will result, unlike the case of v, in 
a decrease of 7 of about one percent for N z > 7. From 
the slopes of these straight lines, we obtain the values of 
7/V. Using the values of v obtained above, we deduce 
the values of 7 which are plotted in Fig. [7| as a function 
of thickness N z . Unlike the case of v, we observe here a 
stronger deviation of 7 from its 2D value (1-75) with in- 
creasing thickness. This finding is somewhat interesting: 
the magnitude of the deviation from the 2D value may 
be different from one critical exponent to another: ~ 3% 
for v and ~ 8% for 7 when N z goes from 1 to 13. We 
will sec below that (3 varies even more strongly. 

At this stage, a natural question arises: does the ab- 



5 




Ln(L) 



FIG. 6: Maximum of susceptibility versus L in the In — In 
scale. 
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FIG. 7: 7 versus N z . 



sence of PBC in the z direction cause these deviations of 
the critical exponents? The answer is no: we have calcu- 
lated v and 7 for N z = 5 in both cases: with and without 
PBC in the z direction. The results show no significant 
difference between the two cases as seen in Figs. [5] and 
[9l We have found the same thing with N z = 11 shown in 
Figs. [TUlandfTTl So, we conclude that the fixed thickness 
will result in the deviation of the critical exponents, not 
from the absence of the PBC. This is somewhat surpris- 
ing since we thought, incorrectly, that the PBC should 
mimic the infinite dimension so that we should obtain 
the 3D behavior when applying the PBC. As will be seen 
below, the 3D behavior is recovered only when the finite 
size scaling is applied in the z direction at the same time 
in the xy plane. To show this, we plot in Figs. 1121 and 
[T31 the results for the 3D case. Even with our modest 
sizes (up to L = N z = 21, since it is not our purpose 
to treat the 3D case here), we obtain v = 0.613 ± 0.005 
and 7 = 1.250 ± 0.005 very close to their 3D best known 
values v ZD = 0.6289 ± 0.0008 and 73D = 1.2390 ± 0.0025 
obtained by using 24 < L < 96) M 

Let us discuss on the deviation of the critical expo- 
nents due to the film finite thickness. For second-order 
transitions, theoretical arguments, such as those from the 
renormalization group, say that the correlation length in 
the direction perpendicular to the film is finite, hence it 



FIG. 8: Maximum of the first derivative of InM versus L in 
the In — In scale for N z = 5 (a) without PBC in z direction 
(b) with PBC in z direction. 
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FIG. 9: Maximum of susceptibility versus L in the In — In 
scale for N z — 5 (a) without PBC in z direction (b) with 
PBC in z direction. The points of these cases cannot be 
distinguished in the figure scale. 



is irrelevant to the criticality, the film should have the 
2D character as long as N z is finite. We have seen above 
that this is not the case here. The deviation begins slowly 
as soon as N z departs from 1. A possible cause for the 
deviation is from the spatially non uniform correlation: 
the correlation in a xy plane depends obviously on its 
position with respect to the surface layer. On and near 
the surface, the spins suffer thermal fluctuations more 
strongly than the interior spins so there is no reason why 
all planes should have the same in-plane correlation be- 
havior even when there is no separate surface transition 
as in the case J s = 1 studied here. Due to this spatially 
non uniform fluctuations, we believe that near the phase 
transition, there are simultaneously several correlation 
lengths which give rise to a kind of "effective" critical 
exponents obtained above. Loosely speaking, we can say 
in another manner that because of its spatial non unifor- 
mity, the correlation in the direction perpendicular to the 
film cannot be separately summed up, it interacts with 
the xy correlation giving rise to " effective" critical expo- 
nents observed in our simulations. In other words, the 
finite thickness makes the dimension of the system some- 
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FIG. 10: Maximum of the first derivative of In M versus L in 
the In — In scale for N z = 11 (a) without PBC in z direction 
(b) with PBC in z direction. 




FIG. 11: Maximum of susceptibility versus L in the In — In 
scale for N z = 11 (a) without PBC in z direction (b) with 
PBC in z direction. 



thing between 2 and 3. Before showing this "effective" 
dimension, we show in Fig. [TJ] the maximum of C™ ax 
versus L for N z = 1,3,5, 13. Note that for each N z 
we had to look for Co, C\ and a/v which give the best 
fit with data of C™ ax . Due to the fact that there are sev- 
eral parameters which can induce a wrong combination 




FIG. 12: Maximum of the first derivative of In M versus L in 
the In — In scale for 3D case. 



FIG. 13: Maximum of susceptibility versus L in the In — In 
scale for 3D case. 
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FIG. 14: ln(C™ ax - Co) versus InL for N z 
The slope gives a/v (see Eq. I12|l . 
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of them, we impose that a should satisfied the condition 
< a < 0.11 where the lower limit of a corresponds to 
the value of 2D case and the upper limit to the 3D case. 
In doing so, we get very good results shown in Fig. [T4l 
From these ratios of a/v we deduce a for each N z . The 
values of a are shown in Table U for several N z . 

It is interesting to show now the effective dimension 
of thin film discussed above. Replacing the values of a 
obtained above in d e Q = (2 — a)/v we obtain (Lg- shown 
in Fig. [[51 

We note that d e Q is very close to 2. It varies from 2 
to ~ 2.061 for N z going from 1 to 13. The 2D character 
is thus dominant even with larger N z . This supports the 
idea that the finite correlation in the z direction, though 
qualitatively causing a deviation, cannot strongly alter 
the 2D critical behavior. This point is interesting be- 
cause, as said earlier, some thermodynamic properties 
may show already their 3D values at a thickness of about 
a few dozens of layers, but not the critical behavior. To 
show an example of this, let us plot in Fig. [T^the tran- 
sition temperature at L = oo for several _/V 2 , using Eq. 
[16] for each given N z . As seen, T c (oo) reaches already 
~ 4.379 at N z = 13 while its value at 3D is 4.51,21 A 
rough extrapolation shows that the 3D values is attained 
for N z ~ 25 while the critical exponents at this thickness 
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FIG. 15: Effective dimension of thin film as a function of 
thickness. 
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FIG. 16: Critical temperature at infinite L as a function of 
the film thickness. 



are far away from the 3D ones. 

We give the precise values of T c (oo) for each thickness. 
For N z = 1, we have T c (oo) = 2.2701 ± 0.0003 from T c 
of specific heat and 2.2697 ± 0.0005 from T c of suscep- 
tibility. From these we have T c (oo) = 2.2699 ± 0.0005. 
Note that the exact value of T c (oo) is 2.26919 by solving 
the equation sinh 2 (2 J/T c ) = 1. Again here, the excellent 
agreement of our result shows the efficiency of the multi- 
ple histogram technique as applied in the present paper. 
The values of T c (oo) for other N z are summarized in Ta- 
ble D 

Calculating now M(L) at these values of T c (oo) and 
using Eq. [T5l we obtain (3/v for each N z . For N z = 1, 
we have 0/v= 0.1268 ±0.0022 which yields (3 = 0.1266± 
0.0049 which is in excellent agreement with the exact re- 
sult (0.125). Note that if we calculate (3 from a+2/3+7 = 
2, then /3 = (2 - 1.75198 - 0.00199)/2 = 0.12302 ± 0.0035 
which is in good agreement with the direct calculation 
within errors. We show in Fig. [T7] the values of (3 ob- 
tained by direct calculation using Eq. [15] Note that the 
deviation of f3 from the 2D value when N z varies from 1 
to 13 represents about 60%. Note that the 3D value of (3 
is 0.3258 ±0.0044,21 

Finally, for convenience, let us summarize our results 
in Table [T] for N z — 1,3, ...,13. Due to the smallness of 
cv, its value is shown with 5 decimals without rounding. 



IV. CONCLUDING REMARKS 

We have considered a simple system, namely the Ising 
model on a simple cubic thin film, in order to clarify the 
point whether or not there is a continuous deviation of 
the 2D exponents with varying film thickness. From re- 
sults obtained by the highly accurate multiple histogram 
technique shown above, we conclude that the critical ex- 
ponents in thin films show a continuous deviation from 
their 2D values as soon as the thickness departs from 1. 
We believe that this deviation stems from deep physical 
mechanisms, not from the calculation method used here. 
We would like moreover to emphasize some additional 
interesting observations: 

1. The deviations of the exponents from their 2D val- 
ues are very different in magnitude: while v and a vary 
very little over the studied range of thickness, 7 and spe- 
cially (3 suffer stronger deviations 

2. With a fixed thickness (> 1), the same critical expo- 
nents are observed, within errors, in simulations with or 
without periodic boundary condition in the z direction 

3. To obtain the 3D behavior, finite size scaling should 
be applied simultaneously in the three directions. If we 
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FIG. 17: Critical exponent (3 versus the film thickness ob- 
tained by using Eq. 1151 



do not apply the scaling in the z direction, we will not 
obtain 3D behavior even with a very large, but fixed, 
thickness and even with periodic boundary condition in 
the z direction 

4. With regard to the critical behavior, thin films be- 
have as systems with effective dimensions between 2 and 
3, depending on the film thickness. Note however that, 
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TABLE I: Critical exponents, effective dimension and critical temperature at infinite xy limit as obtained in this paper. 



N z 


V 


7 


a 


P 


doff 


T c (oo) 


1 
3 
5 
7 
9 
11 
13 


0.9990 ± 0.0028 
0.9922 ± 0.0019 
0.9876 ± 0.0023 
0.9828 ± 0.0024 
0.9780 ± 0.0016 
0.9733 ± 0.0025 
0.9692 ± 0.0026 


1.7520 ± 0.0062 
1.7377 ±0.0035 
1.7230 ± 0.0069 
1.7042 ± 0.0087 
1.6736 ± 0.0084 
1.6354 ± 0.0083 
1.6122 ±0.0102 


0.00199 ± 0.00279 
0.00222 ±0.00192 
0.00222 ± 0.00234 
0.00223 ± 0.00238 
0.00224 ±0.00161 
0.00224 ± 0.00256 
0.00226 ± 0.00268 


0.1266 ±0.0049 
0.1452 ± 0.0040 
0.1639 ±0.0051 
0.1798 ±0.0069 
0.1904 ±0.0071 
0.1995 ±0.0088 
0.2059 ± 0.0092 


2.0000 ± 0.0028 
2.0135 ± 0.0019 
2.0230 ± 0.0023 
2.0328 ± 0.0024 
2.0426 ±0.0016 
2.0526 ± 0.0026 
2.0613 ±0.0027 


2.2699 ± 0.0005 
3.6365 ± 0.0024 
4.0234 ± 0.0028 
4.1939 ±0.0032 
4.2859 ± 0.0022 
4.3418 ± 0.0032 
4.3792 ± 0.0034 



except a strong deviation of 7, other exponents stay near 
their 2D limit even with a large thickness, while non crit- 
ical thermodynamic properties may attain 3D behaviors 
at a thickness of about a few dozens atomic layers. 

To conclude, we hope that the numerical results shown 
in this paper will stimulate more theoretical analysis in 
search for the origin of the continuous variation of the 
critical exponents with changing thickness. It should be 
also desirable to study more cases to clarify the role of 



thickness on the transition behavior of very thin films, 
in particular the effect of the film thickness on the bulk 
first-order transition. 
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